Fidelity of the surface code in the presence of a bosonic bath 
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We study the resilience of the surface code to decoherence caused by the presence of a bosonic bath. 
This approach allows us to go beyond the standard stochastic error model commonly used to quantify 
decoherence and error theshold probabilities in this system. The full quantum mechanical system- 
bath dynamics is computed exactly over one quantum error correction cycle. Since all physical 
qubits interact with the bath, space-time correlations between errors are taken into account. We 
compute the fidelity of the surface code as a function of the quantum error correction time. The 
calculation allows us to map the problem onto an Ising-like statistical spin model with two-body 
interactions and a fictitious temperature which is related to the inverse bath coupling constant. The 
model departs from the usual Ising model in the sense that interactions can be long ranged and can 
involve complex exchange couplings; in addition, the number of allowed configurations is restricted 
by the syndrome extraction. Using analytical estimates and numerical calculations, we argue that, 
in the limit of an infinite number of physical qubits, the spin model sustain a phase transition which 
can be associated to the existence of an error threshold in the surface code. An estimate of the 
transition point is given for the case of nearest-neighbor interactions. 
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I. INTRODUCTION 

Recent progress in implementing controllable multi- 
qubit systems in the laboratory has sparked renewed in- 
terest in topological quantum computing schemes. Par- 
ticular attention has been devolted to the surface code 
1-[3|, which is a planar version of Kitaev's toric code 
J|. From a physical implementation viewpoint, the sur- 
face code has two important advantages in comparison 
to other schemes: (i) all gates are local, and (ii) simula- 
tions indicate that the topological protection yields very 
high tolerance for errors. The latter is based entirely on 
stochastic models for errors. These models point to er- 
ror threshold probabilities per single qubit operation or 
cycle ranging from near 1% [3, Q up to 19% Q- The 
large error threshold comes at the expense of hardware: 
a vast number of local operations and physical qubits is 
required to build a useful computing machine p]. Yet, 
this tradeoff seems attractive nowadays for a number of 
physical realizations such as cold atoms [3] , ion traps T0| , 
Rydberg atoms |ll| . semiconductor systems 12], and su- 
perconducting integrated systems |13| . 

Most error threshold estimates so far have relied on 
the assumption of errors being uncorrelated in time and 
space. However, given the large-scale integration that 
will be required to implement a surface code, this as- 
sumption seems unwarrantable on physics grounds. The 
need to have tens of millions of physical qubits siting on 
a common substrate and interacting with each other and 
with the controlling eletronics is very likely to introduce 
environmental modes, which will effectively couple the 
time evolution of the physical qubits. Under these cir- 
cumstances, errors will become correlated and it is un- 
clear whether the system will retain its high error thresh- 
old. In fact, previous studies of the impact of correlated 
errors on standard (non-topological) quantum computing 
codes have shown that error thresholds may be reduced 



or altogether disappear in some situations [lJ-[3- In- 
vestigating the effect of correlations between errors in 
the surface code is the main goal of this paper. 

In a recent paper [20| , two of us showed that the time 
evolution of the surface code in the presence of a com- 
mon bosonic bath can be mapped onto a statistical spin 
model. This mapping allows for the computation of the 
surface code fidelity much in the same way that one com- 
putes the partition function and expectation values in a 
spin model. As a result, the existence of an error thresh- 
old was related to the existence of a phase transition in 
the statistical model. Even though the interpretation of 
the crossing of the error threshold as a classical phase 
transition is not new [3|, |2l|, our formulation is novel 
since it takes into account the full quantum mechanical 
time evolution of the qubits in the presence of a dynam- 
ical environment. In addition, rather than evaluating er- 
ror probabilities, we compute directly the fidelity of the 
logical qubit. Our choice of environment, a collection of 
freely propagating massless bosonic modes, is realistic for 
systems where decoherence can be related to the coupling 
to phonons, magnons, and electromagnetic modes. 

Below, we provide a detailed derivation of the evolu- 
tion operator of the combined surface code-bosonic bath 
system. We focus our attention on a single quantum er- 
ror correction cycle and assume that, after the syndrome 
extraction, the bath is reset to its ground state. Within 
this approximation, we find that the fidelity can be writ- 
ten as a function of the expectation value of single-qubit 
logical operators. The study of these expectation values 
can be related to the physics of an Ising-like spin model 
with a complex fictitious temperature. We use both exact 
and mean-field finite numerical calculations to argue that 
the spin model sustains a thermodynamic phase transi- 
tion in the limit of an infinite number of physical qubits. 
The critical temperature of the spin model yields a cou- 
pling constant threshold value, which is found to depend 



mainly on bath parameters. 

The paper is organized as fohows. In Sec. |lT]we give a 
give a brief introduction to the essential elements of the 
surface code and set some of the notation used later. Sec- 
tion mil presents a Hamiltonian formulation of the prob- 
lem in terms of bosonic modes coupled to physical qubits 
which allows us to obtain a compact form for the evolu- 
tion operator of the combined logical qubit-bath system. 
The evolution operator involves a bath correlation func- 
tion which is explicitly evaluated for three representative 
situations. The effect of syndrome extraction on the evo- 
lution operator is described in Sec. llVl and an expression 
for the fidelity in terms of expectation values involving 
qubit operators is derived in Sec. El The mapping of the 
fidelity calculation onto a statistical model is given in 
Sec. I VII and the connection between the fictitious critical 
temperature and the error threshold probability is shown 
in Sec. IVIII In Sec. IVIIII we estimate the fictitious critical 
temperature via a low-temperature expansion. Numer- 
ical supporting the existence of an error threshold are 
shown in Sec. IIXI which is a very encouraging result. 
Conclusions and a critical discussion of the approxima- 
tions involved and future directions of investigation are 
drawn in Sec. [Xj A number of appendices with technical 
details of the calculations are also provided. 



II. SURFACE CODE 

Following Ref. HI, we define the surface code as collec- 
tion of N spins 1/2 (physical qubit systems) located on 
the edges of a two-dimensional lattice with two types of 
boundaries, as shown in Fig. [T] The lattice comprises n 
and m qubit rows and columns, respectively. Measure- 
ments are done on two types of stabilizer operators: stars 
A(^, which are associated to lattice vertices (0), 



Ac^ 






(1) 



and plaquettes B\j, which are associated to tiles (D), 
including the ones at the open boundaries. 



Br 



n 



cr, . 



(2) 



In Eqs. (HI) and ©, the Pauli spin operators (Ji act on 
qubits. Thus, there are A^n = {n + l)m plaquette opera- 
tors and A'^o = {m + l)n star operators. The A^ physical 
qubits store one logical qubit. There are nz, = 2 distinct 
logical operators: X and Z. They are formed by a string 
of physical qubit operators along paths that cut through 
the lattice: 



and 



ieTz 



(3) 



(4) 



where Tz runs between qubits at opposite open bound- 
aries (left to right), passing through vertices along the 
way, while Tx runs between qubits at opposite closed 
boundary (top to bottom), crossing tiles (see Fig. [T|). 
Notice that vertices and tiles form dual lattices. 
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FIG. 1. (Color online) A 3 x 3 two-dimensional square lat- 
tice structure for a surface code. Physical qubits (circles) are 
located at the edges of the lattice, which has open (vertical) 
and closed (horizontal) boundaries. In the general case, the 
lattice is n X m in size, with nm vertical and (n -I- l)(m + 1) 
horizontal edges and N = physical qubits. The light colored 
strips show possible paths for the logical operators Z and X. 
a is the lattice constant. 

The protected code space contains two states, |t) and 
|J). Both states are eigenstates of all stabilizer operators 
with eigenvalue -1-1. Errors can be inferred by measuring 
the stabilizer operators and tracking down which stars 
or plaquettes yielded —1 values. A decoding procedure 
is needed to decide which recovering operation to perform 

M- 

The code state can be generated by the action of a 
product involving all star operators on the z ferromagnet 
state, namely. 



and 



where 



and 



G = 




N 



m = l[\t)^.. 



(5) 



(6) 



(7) 



(8) 



ler^ 



Notice that the product in Eq. ([7]) can be expanded as 

n(i+^o)=i+E^^+ E ^oi^o.+---- (9) 

0i#02 

The number 2^* appearing in the prefactor of G is the 
number of terms appearing in the expansion: 2^* = 

ro*)+(t)+---+(;^:)- 



III. THE BOSONIC ENVIRONMENT 

The most general bath model would allow for both flip 
and phase errors to occur. However, only perturbative 
calculations would be possible in this general case. Since 
our goal is to obtain nonperturbative results, we focus our 
discussion on flip errors only (it is possible to rephrase the 
model to induce only pure dephasing model by a simple 
change of basis). The Hamiltonian we consider is written 
as 



H = H„ + V, 



where Hn is a free bosonic Hamiltonian, 



Hn 



22 "^k a^ak, 

k/O 



and 



V 



2 ^^ 



fin 



(10) 



(11) 
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where r^ denotes the spatial location of a qubit i and / 
is a local bosonic operator. 



/(r) = 



{v/uq) 
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flk 



(13) 

Here, D is the bath spatial dimension, v is the bosonic 
mode velocity, Wk = w|k|, and luq is a characteristic fre- 
quency of the bath (notice that / is dimensionless since 
we adopt units such that H — 1). The creation and anni- 
hilation operators of the bosonic modes follow the stan- 
dard commutation relations, namely, [ak,aj^,] = (5k, k' and 

K>«k'] = [aL«k'] =0- 

The choice of s depends on the physical nature of the 
environment and on which bosonic degree of freedom cou- 
ples to the qubits. When the qubits couple directly to the 
bosonic displacement field, we choose s = —1/2, whereas 

I 



when they couple to the bosonic current operator, we 
choose s = 1/2 instead. Notice that these two choices 
allow us to write the coupling between the qubits and 
the bosonic environment as a simple function of the free 
bosonic field. Hence, they render an interaction Hamil- 
tonian with commutators that are subluminal, namely, 
which are equal to zero outside the boson light cone (see 
below). But this is not the general rule. For instance, 
a model that couples the environment to two different 
qubit components would render an interaction Hamilto- 
nian with non-subluminal commutators (regardless of our 
choice for s). 

This apparent problem comes from the fact that we 
usually think of errors in a dynamical sense: they happen 
in a point in space-time and create bosons that propa- 
gate at the speed of light. However, this is an incorrect 
interpretation to the equations we are about to derive. 
We will not be considering pulses propagating thought 
a medium, but rather looking at allowed normal modes 
of the bath and how they relate to different qubits con- 
figurations. If there is no fundamental symmetry rea- 
son for their suppression, long- wave length modes of the 
bath will in general introduce superluminal effective in- 
teractions. A very simple way to highlight this fact is 
to rewrite the bosonic model in a coherent state basis, 
flk = Ok + ck, where au and aj^ also obey standard com- 
mutation relations and a in a constant. This procedure 
introduces an effective instantaneous interaction between 
the qubits as much as the Coulomb gauge introduces the 
instantenous Coulomb interaction in quantum eletrody- 
namics. 

For a time interval A, the error model comprised by 
Eqs. PH)) to (|12p leads to the following evolution operator 
in the interaction picture [l6| : 



C/(A) = Tt exp 



.A 



.A 

/ dtY,fir^,t)a^ 
Jo 



(14) 



Combining a Magnus expansion with the Zassenhaus for- 
mula (see Appendix |X]), we arrive at a remarkably simple 
expression for the evolution operator. 



C/(A) = X cxp 



-TE*-,(A)afa| 



»#j 



: exp 



iX 



E^r.(A)af 



(15) 



r 



where 



cI>rs(A) = -[t;(f)(A) + e(,^)(A) 



X = exp 



E*-r.(A) 



(16) 

(17) 



and :: stands for normal ordering. In Eqs. (J15p and (1161) . 



we have introduced two bath correlation functions, 
g«(A) = i / dh r dt2{[J{vM)J{^M)] 



and 



+ [/(s,ii),/(r,t2)]} 



t;(f)(A) = (0|Fr(A)F,(A)|0), 



(18) 
(19) 



and the auxiliary bosonic field 



Fr(A) 



A 



dtf{r,t). 



(20) 



(Notice that x is a real number since Qrr — 0.) 



where 0{x) denotes the Heaviside step function. The 
imaginary part reads (see Appendix IC 1|) 



G^iH^) 



(uA- |r-s|) 




(22) 



Below, we present the functional form of the correla- 
tions functions for two-dimensional bosonic baths {D — 
2). Details of the calculations are provided in Appendices 
IB] and \C\ We consider three representative values of the 
power s which appears in the qubit-bath coupling con- 
stant dependence on the bath mode momentum [see Eq. 
(IT^ ]. These values are s = —1/2, 0, and s = 1/2, cor- 
responding to subohmic, ohmic, and superohmic baths, 
respectively. This classification is standard and follows 
from the bath's spectral function frequency dependence 
at low frequencies: sublinear (subohmic), linear (ohmic), 
and superlinear (superohmic) '22|. 



Subohmic bath 



For two-dimensional subohmic baths (D = 2 and s = 
— 1/2), the bath correlation function takes a simple closed 
form. The real part reads (see Appendix IB ip 
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(21) 



B. Ohmic bath 



Choosing D — 2 and 5 = 0, the real and imaginary 
parts of the bath correlation function take the forms [23| 
(see Appendices IB 21 and IC2|) 



gi^H^) 



■ arcosh 



vA 



and 



sii'i^)-^ 



TTiOn 



-e(vA 

2 ^ 



(vA 



s|) 



s|), (23) 



vA 



e{\r 



vA) 



.(24) 



Notice that the real part of the correlation function 
vanishes for distances larger than vA. For this bath as 
well as others, the number of qubits within the spatial 
range of the correlation function is determined by the 
ratio f A/a, where a is the surface code lattice constant. 



C. Superohmic bath 

Choosing D — 2 and s = 1/2, we find for the real part 
(see Appendix IB 3P 



Gi^\A)^ 



9{\r- 



vA) 



r- s 



V^|r - s|2 - (t;A)2 

For the imaginary part we find (see Appendix IC 31 

iv ^(vA-|r-s|) 



GiiH^) = 



7ra;3 ^(wA)2-|r-s|2' 



(25) 



(26) 



A schematic representation of the spatial dependence 
of these correlations functions is shown in Fig. [2] 



IV. SYNDROME EXTRACTION 



Let us assume that the system is prepared initially in 
the logical state If) and the boson field initial state is 



the vacuum, 



\^o) = {Gm)^\o)b. 



(27) 



We then let the system evolve under the interaction 
Hamiltonian until a time A, when an error correction pro- 



\s = -m\ G('^> 




FIG. 2. Schematic representation of the spatial dependence of the correlation functions CJis (A) and tj^s (A) for s 
and 1/2. We use d= |r — s|. 
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tocol is performed flawlessly. The syndrome extraction 
operator is equivalent to the application of the projector 



P 



1 



2^0 2^0 



l[il + baBa) II il + aoA(^), (28) 



where a\j — ±1 and &^ = ±1 are the syndromes for each 
star and plaquette operator, respectively. Since we are 
assuming only bit-flip errors, the projection over stars is 
just the identity for ao = 1, namely, 



^=in(l + ^n^n) 



(29) 



Using that ^o^ = C! and [X, G] — 0, we can rewrite the 
projector in a slightly different form. 



P^R\t){t\R + R\l){l\R, 



(30) 



where R is the recovery operation chosen to be per- 
formed. 

In principle, we should consider all possible syndrome 
outcomes. However, it is useful to look at the most be- 
nign evolution and assume 6n = 1. Such an evolution 
corresponds to the system remaining in the vacuum of 
the gauge fields after a time A. This nonerror syndrome 
provides an upper bound to the available computational 
time. In addition, it simplifies the calculations by remov- 
ing from consideration a recovery operation that tries to 
steer the system back to the computational basis [2J, [2^ . 
Thus, for this particular case. 



p = |t)(t| + |;)(;| 



(31) 



The environment state is unaffected by the error cor- 
rection protocol. If no extra step is taken to dissipate 



excitations that pile up over time, the environment will 
keep a memory of events that happened during the QEC 
period. Keeping track of such excitations between QEC 
cycles in a fidelity calculation is a difficult task even for 
simple, non-topological logical qubit systems [la. Im [l8| . 
For topological qubits, the task is considerably harder 
due to the exponential number of terms that enter in the 
composition of the computational states. 

Thus, in order to make the formulation amenable to 
an analytical calculation, we consider an extra step to 
the QEC protocol. In addition to projecting the quan- 
tum computer wave function back to the logical Hilbert 
space, we assume that at the end of the QEC step the 
environment is reset to its ground state. This is equiva- 
lent to imposing at the end of the QEC step the projector 
limj'b^tij^o 6"^°''^^"^'"'"", for some environment tempera- 
ture Tbath defined with respect to some even larger reser- 
voir. 

A consequence of this extra QEC hypothesis is that 
we exclude from the calculation any spatial correlation 
between QEC periods, as well as memory and spatial 
correlations between the time evolution of bras and kets. 
This new projector operator can be conveniently written 
as 



P' = |*o)(*o|+^|*o)(*o|^. 



(32) 



After the projection, the wave function must be normal- 
ized again. For this purpose, consider the normalization 
factor 



where 



(vI'o|f/t(A)P'L/(A)|*o) = |^r 



^=(*o|C/(A)|*o) 



m 



(33) 



(34) 



and 



and 



B={^o\XU{A)\^o). 



(35) 



Below, we use the expectation values A and B to compute 
the surface code fidelity after one QEC cycle. 



V. FIDELITY 

The fidelity of the surface code after one QEC cycle 
can be defined as 



F= |(*QEC|*0) 



(36) 



where j^I^o) is the initial state of the qubit system and 
the bath and 



|^QEc)=-P'C/(A)|*o) 



(37) 



The expectation values A and B now come in handy since 
they allow one to obtain a simple expression for the fi- 
delity, 



F = 



(^0 |P'^(A)| ^o) (^0 \W{^)P'\ ^o) 
(*o|f/t(A)P'C/(A)|*o) 

(*o|C/(A)|*o)(*o|f/^(A)|*o) 



1/2 



\A\' + \Bf 



1/2 



1^1 



\A\' + \B\' 



1 + 



Mi 

l-A|2 



(38) 



where we used that P'\^q) = |*o)- Thus, our task of 
determining the fidelity is reduced to evaluating the ratio 

\B\y\A\\ 



/3 = 



27r 



1 



(^oA) 



_D+2s-2 ■ 



(42) 



Clearly, Eq. (|4T|) can be interpreted as an effective in- 
teraction between qubits intermediated by the environ- 
mental bosons. The connection to a statistical model 
becomes more apparent when we consider that the ferro- 
magnetic state along the z direction can be expanded in 
in the x basis, namely. 



N 

n 



I t)^,x + I i)^.a 

V2 



(43) 



Inserting this expression into Eqs. ([M| and (P(7|. we 
arrive at 



^=WT.^'''''(S\G'\S) 



and 



s 



(44) 



(45) 



where S stands for the eigenstates of the operator 
^^=l < and 



(3Es = {Sm\S). 



(46) 



Notice that the expectation values in Eqs. (H^ and pS)) 
vanish for those states \S) where at least one start oper- 
ator has a —1 eigenvalue. Therefore, those equations can 
be rewritten as 



■^~ ON Z^ 



'PEs, 



(47) 



VI. MAPPING ONTO A STATISTICAL MODEL 

Let us first consider the cases where the bath correla- 
tion function $rs(A) has both real and imaginary parts 
finite, namely, < D -|- 2s < 3. The insertion of the 
evolution operator given by Eq. ([T5|) into Eqs. ([M|) and 
([55]) results in very compact expressions for A and B 
since these expectation values are taken on the bosonic 
vacuum. After a short manipulation, we arrive at 



and 



A = x{Fz 



-pn 



G^\F,) 



and 



where we introduced 



(39) 



(40) 



(41) 



^ S' 



■I3E, 



■' {s'\x\s'). 



(48) 



where the sums are over the subset of states S' where all 
star operators take positive values: 



{S'\Ao\S')^+l. 



(49) 



It is clear now that A is proportional to the partition 
function of a classical statistical spin model with a re- 
stricted configuration space. Then, C is equal to the ex- 
pectation value of the operator X in this model. 

The statistical model defined by Eqs. (gl]), (gS]), (|47| . 
and (|48p is nontrivial in a number of ways. First, the 
interaction term (|4ip is not purely real. Second, the in- 
teraction range is not necessarily restricted to nearest 
neighbors. Third, the constraint imposed by Eq. (j49p 
severely reduces the size of the configuration space. 



Since (S" \x\ S') can only take the values ±1, one can 
rewrite Eqs. (|47l) and pS)) as sums over "energy" eigen- 
values, namely, 



^=^EL9+(^')+5-(i?')] e 



-/3E' 



(50) 



strength of the coupling between the qubits and the envi- 
ronment. The effective exchange interaction amplitude, 
Jij, depends on the bath characteristics (e.g., spatial di- 
mension and spectral density), on the QEC cycle dura- 
tion A, and on the ratio a/(Aw). 

For instance, consider the ohmic bath, where 



and 



E' 



-I3E' 



(51) 



where g^{E') are the number of qubit configurations with 
energy E' and {X) = ±1. The prime indicates that only 
configurations where all star operators have -1-1 expecta- 
tion value are considered, Eq. (|T7)) . 

When the sums in Eqs. ([50]) and (|51l) are not restricted 
by Eq. (|47p . the time-reversal symmetry of the Hamilto- 
nian implies g~^{E) = g^{E) for logical operators X con- 
taining an odd number of af qubit operators. Therefore, 
in this case, B = 0. For X containing an even number 
of af operators, for each "energy" eigenvalue E, {X) is 
either +1 [and g-{E) = 0] or -1 [and g+{E) = 0], but 
the value of B cannot be predicted. 

In the case of a restricted sum, it is straightforward to 
see that the separation of configurations in time-reversal 
symmetry classes is not useful. Consider that at the ver- 
tical boundaries one can form operators A() with three 
qubits. In this case, even if a certain configuration |S") 
satisfies Eq. (|47p . its time- reversal partner will not and 
therefore will not be included in Eqs. ([50]) and ([5T|). 
Thus, the restriction is equivalent to projecting out time- 
reversal partner of |5"). 

As explained in Ref. 120, one useful way to understand 
this point is to break up the states {|S")} into two groups, 
(IS*;)} and {\S'_}}, where 



i^;) = II^ai^- 



(52) 



/3 



1 

2^ 



(54) 



Using Eqs. (p3)) and ([24)) . the effective Hamiltonian of 
the statistical model can be written as 



with 



n 



arcosh 



•^u' ^ 2 ^ 



I arcsm 






\\ri-r,\J 



\r.- 



A 



vA 
vA 



<1, 
> 1. 



(55) 



(56) 



The real part of Jij is nonzero only between qubits 
within the light cone of the bosonic modes. The imagi- 
nary part is present for any pairs of qubits, but decays 
rapidly (approximately with the cube of the inverse dis- 
tance) when qubits are outside the light cone. For a lat- 
tice of size L, an extreme limit occurs when vA ~ L, in 
which case all qubits are correlated. In the opposite limit, 
when the QEC cycle period A is sufficiently short (or, 
equivalently, that the lattice constant is large enough), 
so that a/\/2 < vA < a, only qubits belonging to the 
same plaquette are within the light cone defined by the 
free spatial propagation of the bosonic modes and the 
real part of the coefficient J^ vanishes beyond nearest 
neighbors. 



VII. CONNECTING /3 TO THE QUBIT ERROR 
PROBABILITY 



and 



\S') - Zr\S' 



(53) 



Here, Y[j Bo is a product of all plaquette operators that 
do not touch the logical error Zr for a given path F. It 
is then possible to show that this separation leads to the 
appearance of an effective local magnetic field that acts 
only on the qubits along the path F. This local magnetic 
field leads to the time-reversal symmetry breaking in the 
computation of the expectation values in Eqs. ((50)) and 

(EH). 



A. Effective interaction and fictitious temperature 

The parameter /3 plays the role of inverse temperature 
in the statistical model. From Eq. (H^ . we see that 
/? is proportional to A^ , thus serving as a measure of the 



It is useful to relate the fictitious inverse temperature 
/3 of the spin model to the probability p of a qubit hipping 
its spin state during the QEC cycle. The latter can be 
defined as 



P 



(t,|C/|(A)| ;,)(;, |[/,(A)|t,)®|0), (57) 



where {| ti)i I ij)} a-re states of the qubit located at r^, 
|0) is the bath ground state, and 



C/j(A) = rtexp 



-^^ di/(r„t)aj 



(58) 



is the evolution operator of that qubit coupled to the bath 
when the dynamics of all other qubits is frozen. The steps 
in the evolution of Eq. [57] are similar to those used in 
the derivation of the fidelity. The details are provided in 
Appendix iDl The result is 



^' 



exp 



'T^^^(^) 



(59) 
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Notice that for X — 0, p = 0. As the coupling be- 
tween qubits and bath grows in magnitude, p approaches 
1/2, which signals a complete randomization of the qubit 
state. 

The functional relation between p and the fictitious 
temperature (3 can be easily established by evoking Eq. 
(|42|) and employing the explicit form of t/^^r^ (A) as given 
in Ec^. (JB1[) . One obtains 



ln(l-2p) = - 



tt(3{vA) 



D+2S-2 



LD 






-[l-cos(|k|i;A)]. 

(60) 

The sum over momentum diverges in the ultraviolet when 
D + 2s > 2 and the relation between the error probabil- 
ity and the fictitious inverse temperature becomes cutoff 
dependent. For instance, for D = 2 and s — (ohmic 
bath), one finds P= \[^- (2uAA)~''/^] , where A is the 
ultraviolet momentum cutoff. However, for D — 2 and 
s = —1/2 (subohmic case), one finds P = | (l — e"^^/'*), 
which is cutoff independent. 



VIII. ESTIMATE OF /3c 

We now provide an estimate of the critical inverse fic- 
titious temperature, takiiig a slightly different approach 
from that used in Ref. [20. Let us consider the case 
when the effective spin coupling in Eq. (j4ip is real and 
only nearest-neighbor interactions occur, $rir (A) = J 
for |ri — Tj I < a/\/2 and $rirj (A) = otherwise. We will 
carry out a low-temperature (large /3) expansion of Eqs. 
([SH]) and ([5T|) . The key element we exploit in this expan- 
sion is the following property of the surface code: At the 
boundaries of the surface code, the star operators are de- 
fined by three qubits instead of four. This means that if a 
certain state \S) is included in the restricted sums defin- 
ing A and B, its time-reversed counterpart is included. 
This is because reversing all the qubit of a configuration 
with an eigenvalue -1-1 for all star operators yields a —1 
eigenvalue for the star operators at the boundaries. 

This property is particularly useful in the limit of /3 — >■ 
c» when the term with the minimum energy, e"^^™™, 
carries the leading contribution to the sums. For /3 — ?• oo 
and for J < 0, the only state with minimum energy is 
a ferromagnetic state where all spins are pointing along 
the positive x direction, \Fx). Since the the ferromagnet 
state with spins pointing along the negative x is not part 
of the restricted sum, g~{E'^^^ = 0. In this limit, the 
spin model is in the ordered phase, with A = B^ resulting 
in J" = 1/2 (lost fidelity). 

Consider now a large but finite /?. Starting from the 
state l-Fx), the states appearing in A and B can be sep- 
arated into two groups, as shown in Eqs. ([5^ and ([5^. 
The first group, {5'^}, corresponds to the states counted 
in <;+(£"), whereas the second group, {-S*^}, is accounted 
for by g~{E'). Therefore, g~{E') is equal to the number 
of states of energy £" with a logical error Zy for a given 



path r. At large /3, the energy cost of these states is of 
the order of the length of Zr and they are suppressed in 
comparison to other states. The leading terms contribut- 
ing to the sums are the minimum energy state \Fx) and 
the states \S'^) containing only small loops. However, as 
the system size increases the multiplicity factor g^{E') 
increases as well. Its value is proportional to the number 
of ways one can apply the Z operation, or equivalently, to 
the number of self-avoiding walks (SAW) from one open 
boundary to its opposite. The number of SAWs with a 
length I, is related to connective constant /i of the lattice 
and scales as /i'. 

If the multiplicity factor g~{E') is high enough, it can 
compensate the Boltzmann factor suppression. Then, as 
/3 decreases to a certain value /3c, for some energy E'_^ 
a term of the type ^'z(^*) e~^''^* will appear in B with 
the same order as the leading term related to the \S'_^) 

states, namely, e"^"^^-"". Here, 1 2 is the length of Z. 
This criterion provides a crude estimate for Pc'- 



leading to 



-P-^E'^ 



lie 



M'2(^*)e- 



P'=K 



EL - EL„ ' 



(61) 



(62) 



The difference between i?^ and E'^^^ is proportional to 
the length of the logical error Z. Then, the denominator 
is of the order of 2n IzJ, where n is the number of qubits 
interacting with the qubits comprising the logical error Z 
through the Hamiltonian in Eq. ([55]) . A range of possible 
values for the connective constant of a square lattice can 
be found in the literature. If we adopt n = 2.64 j2^], set 
n = 4, and insert these values into Eq. (|62)). we obtain 
/3cJw0.12. 



IX. NUMERICAL EVALUATION OF THE 
FIDELITY 

In light of Sec. |Vll we can interpret the ratio B/A as 
the expectation value of the Mx = Sien '''J' namely, the 
X magnetization of a linear set of spin 1/2 particles em- 
bedded into a spin system governed by the Hamiltonian 
of Eq. ([55]) with the restriction imposed by Eq. (|47p . 
(Here, VI denotes the path defining the logical error X.) 
In the absence of such a restriction, the computation of 
B/A in the thermodynamic limit would follow standard 
procedures used in Statistical Mechanics. The restric- 
tion, however, makes an analytical computation rather 
difficult, if not impossible. Therefore, we resort to nu- 
merical calculations, both exact and approximate, to find 
how B and A (and thus the fidelity F) behave as a func- 
tion of /3 and how this behavior scales with increasing 
system sizes. 

Below, we focus on the case where the effective in- 
teraction strength Jij is real and only involves nearest- 
neighbor qubits. As mentioned earlier, this special case 
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FIG. 3. (Color online) Surface code fidelity of code spaces of 
25 and 41 physical qubits in contact with a bosonic bath when 
star operators are restricted to positive values (^o = 1). 



is of significance to experiments where vA is of the order 
of a. Short-range correlated errors in this case can be in- 
troduced by any measurement or operation on individual 
stars and plaquettes. 



A. Exact calculations 

For two lattice sizes, A^ = 25 and 41, we computed 
A and B for an X operator that ran vertically through 
the middle of the lattice. The computation was done by 
exhaustive enumeration of all orthogonal qubit configura- 
tions \S) that complied with the constraint {S\G^\S) ^ 0, 
namely, that produced only positive plaquette eigenval- 
ues. We verified that the results were insensitive to the 
choice of operator X . The resulting fidelity is shown in 
Fig. [3] as a function of the inverse fictitious temperature 
/3. For small /3 (equivalent to small coupling constant 
A), the fidelity stays close to 1 after one QEC cycle. As 
/3 increases, the fidelity decays and tends asymptotically 
to l/sqrt2, which is expected when B = 0. Another 
important feature is that the transition from J-" ~ 1 to 
J- = l/\/2 becomes sharper as the system size is in- 
creased. This is the expected behavior when, in the ther- 
modynamic, infinite-size limit, a phase transition occurs 
at some critical value of /3. 

We have tested that this behavior is not substantially 
altered when the coupling constant J gains a constant 
imaginary part. The results are shown in Fig. 21 The 
main effect of adding an imaginary part is to create os- 
cillations in the decay of the fidelity as a function of /3. 
The larger the magnitude of the imaginary part in J, 
the more oscillations are observed. However, the relative 
amplitude of these oscillations decrease with incresing 
system size. In the limit of a large number of physical 
qubits, we expect that the oscillations to be relatively 
small and concentrated near the critical value Pc- 

In order to determine the critical value /3c, we resort 
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FIG. 4. (Color online) Fidelity of a code space of 25 physical 
qubits in contact with a bosonic bath when star operators are 
restricted to positive values {A(^ — 1) and an imaginary part 
is added to the coupling constant: J — Jr + iJi. The data 
sets correspond to different values of Ji. 



to the coherent anomaly method, which has been exten- 
sively and successfully used to determine critical temper- 
atures in interacting spin systems [23, [23] . 



B. Mean- field solution: coherent anomaly method 

In the coherent anomaly method (CAM), a cluster 
of interacting spins is embedded inside a mean-field 
medium. Self-consistency is obtained by allowing the 
spins at the boundary of the cluster to experience the 
mean field, which is set equal to the mean value of the 
central spin in the cluster. This constraint provide an 
equation from which the critical temperature can be de- 
termined. As the cluster size is increased, the expecta- 
tion is that the critical temperature obtained in this way 
rapidly approaches the exact value of an infinite-size sys- 
tem [23, [ii]. 

More precisely, let ^g denote the central spin operator 
and let "H be a Hamiltonian describing nearest-neighbor 
interactions inside the cluster, 'Hci, as well as the action 
of an effective field 4>cs at the boundary spins, 



iedn 



(63) 



where dQ denotes the cluster boundary. The expectation 
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value of the central spin is given by 



5.5 



(^o) = 



Tr [^0- 



-/3W1 



Tr [e-^«] 



(64) 



where the trace is carried over all allowed spin configu- 
rations. Expanding the exponentials in Eq. ([63|) to the 
lowest nontrivial order in the effective field, we find that 



(5o^) = (5i?)ei-/3J0off V(5o"Si)ci 



iedn 



(65) 



where (■■■)ci denotes the expectation value taken with 
just the Hamiltonian "Hci and neglecting the boundary 
field. Setting (Sq) equal to 0off , we arrive at the equation 



i-i3,jj2{sssr. 



0, 



(66) 



iedn 



which can be solved numerically to yield the critical in- 
verse temperature /3 = /3c as a function of J. The most 
costly part of the procedure is the calculation of the cor- 
relation function {S§Sf)ci, which requires an exhaustive 
enumeration of all spin configurations within the cluster. 
We employed this method to compute the critical value 
of /3 for surface codes with clusters of increasing sizes and 
performed a finite-scaling analysis to estimate the critical 
value in the thermodynamic limit. The result is shown 
in Fig. [5j As in the case of the exact numerical cal- 
culations, we employed the constrained nearest-neighbor 
Ising model of Sec. I VII with J real and only allowed 
for spin configurations with positive plaquette eigenval- 
ues. We find that (^cJ = 0.193(2) for an infinitely large 
system, which is about 60% higher than the estimate 
presented in Sec. I Villi (given the roughness of the ap- 
proximations involved in the estimate, the discrepancy 
seems quite acceptable). The extrapolated value also 
matches quite closely the point where the downturn of 
the fidelity develops (see Figs. [3]andHl), providing addi- 
tional support for the existence of a phase transition in 
the thermodynamic limit in the case of nearest-neighbor 
interactions. 



X. DISCUSSION AND CONCLUSIONS 

We have presented a fully quantum mechanical calcula- 
tion of the fidelity of the surface code in the presence of a 
bosonic bath. We considered the fidelity after a complete 
quantum error correction cycle and in the most benign 
case, when a nonerror syndrome occurs. A important as- 
sumption made in the calculation was the resetting of the 
bath to its ground state after the syndrome extraction. 
We then expressed the fidelity as a function of the ra- 
tio of two complex amplitudes which were formulated as 
expectation values of a statistical spin model with com- 
plex two-body interactions and a restricted configuration 
space. We presented both analytical estimates and nu- 
merical evidence that the statistical spin model sustain 



T //= 5.17(6)- 2.46(22) /L 
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FIG. 5. Finite-size scaling of the critical fictitious tempera- 
ture Tc obtained from cluster mean-field calculations for lat- 
tice of sizes 13, 25, and 41. A real Ising interaction of strength 
J involving only nearest neighbors was used. The circles are 
the numerical data and the dashed line is a linear fit. 



a ordered/disordered phase transition in the limit of an 
infinite number of qubits. The existence of such a phase 
transition can be directly related to the existence of a 
threshold on resilience of the surface code: provided that 
the bath coupling constant remains lower than a critical 
value, fidelity can be maintained close to unity with in- 
creasing system size. This is very good news for those 
interested in large-scale implementations of the surface 
code. 

This work provides more detailed derivations and 
deeper analyses than Ref. l20l in additional to new nu- 
merical support to the existence of an error threshold in 
the surface code in the presence of correlated noise. We 
note that our approach differs substantially from the sur- 
face code literature since we do not rely on a stochastic 
error model. 

In order to understand the difference between a thresh- 
old due to correlations and the usual stochastic model 
discussion, let us consider the case of D = 2 and s = 
(ohmic bath). As shown in Sec. I VIII the probability 
of such a bosonic bath to produce a bit flip error is 
p = i [l - (2wAA)-'3/2j jf ^g tj^j^g ^jje ultraviolet cutoff 
to infinity, then we are bound to find p = 1/2. However, 
in most physics systems, and condensed matter systems 
in particular, the cutoff is finite. Hence, we can expand 
p for small /3 and A to obtain 



/3^ 



4p 



ln|2i;AA| 



(67) 



For the ohmic model, we also found that the real part of 
the effective interaction is given by [see Eq. (f56| ] 



<^7,^ 



im 

2 



vA 



(68) 



for qubits within the causality cone. In both equations, 
the logarithms are slow growing functions and should be 
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regarded as producing numbers of the same order. There- 
fore, we can rewrite p as 



2p 
J' 



(69) 



where we took J ^ Jij ^ (1/2) In |2i;AA|. Now, we have 
also found that the inverse critical temperature is [see 
Eq. dSH)] 



^^ 2nJ' 



(70) 



To be resilient to correlated errors, we must require the 
system to be above the critical temperature, thus /3 < /?c- 
Using the equations above, we find that 



P< 



In/x 
4n ' 



(71) 



For the least correlated case, where only nearest neigh- 
bors effective interactions take place, we have 



P< 



In/i 



6%. 



(72) 



That is, the threshold for the surface code in the presence 
of an ohmic bath is reduced to at most pc ~ 6% due to 
the introduction of nearest-neighbors correlated errors. If 
we allow for longer-range correlated errors, the threshold 
will steadily decline. 

Within our formulation, exact analytical calculations 
of the threshold based on the fidelity are daunting. Thus, 
it is likely that quantitative results will always require 
numerical simulations. We are in the process of simulat- 
ing statistical spin models with more general interactions 
than the nearest-neighbor case investigated here. In ad- 
dition, further investigations are necessary to relax the 
assumption of bath resetting and to evaluate the effects 
of residual qubit correlations on the fidelity over multiple 
cycles. Such studies are also under way. 
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Appendix A: Evolution operator 

Consider Eq. ((T^ . where the bosonic field in the in- 
teraction picture reads 



/(r,i) 



(v/uio) 



D/2+S 



L^/2 



Ei^r 



s / ik-r-j|k|i.t t 



k^^O 



^-^k■r+^\k\vt^^ 



(Al) 



We can write U{A) — exp[J7(A)], where ri(A) follows the 
Magnun expansion 



17(A) = f]i(A) + rJ2(A) + nsiA) 



with 



r!i(A) 



n2(A) 



iX 
'^ Jo 



1 f\ 



dtY^f{v,,t)al 



(A2) 



(A3) 



2 „A 



dti I dt2 
Jo 



2! V2, 

5^[/(r„ti),/(r,,t2)]afaj=, 



(A4) 



hO 



xi;(l/(ri,fi),[/(r,.i2),/(r,,ij)|] 

+ [/(rfe,t3),[/(r„i2),/(r„ii)]]) 

X <(7j(T^, (A5) 

etc. Since [/(r^, ti), /(r^, ^2)] is a c-number, only the first 
two terms in the expansion survive and we can write 



C/(A) = cxp 



A 



-z-^Fr,(A)ar 



X exp 



-tE^^:^.(a)<-i 



2J 



(A6) 



where Qrilj^A) and -Fri(A) are defined in Eqs. (fTS]) and 
(P0| . respectively. 

It is convenient to rewrite the first exponential in Eq. 
(jA6p as a normal ordered term. For this purpose, we note 
that Fr- (A) has the form 



■Pi-. (A) ^^{g*r^Mav, + gr^.yai)^ 



(A7) 



where 



9r,k 



_■iv^^^^_^^^^s^l^^k■r(^^^\k\vA_^y (A8) 
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Thus, we can write 



iX 



exp 



where 



E^r.(AK 



J]^ exp(a|{^Wk-akWk), (A9) 



iX 



Vk 



Ilffr 



,k CT,; 



(AlO) 



Now we can use the Zassenhaus formula appUed to 
bosonic operators, 



This is integral is convergent provided that < D + 2s < 
4. Assuming D = 2, we can write 



5i^\^)- 



2(wM) 



2+2s 






|k|2s-2gik-(r-s) 



X [l-cos(|k|z;A)] 



\2+-2s (-A 

-^^5 / dkk^'-^Jo{k\r-s\) 

71" W^ Jo 



X [1 - COS (kvA)] , (B2) 

where Jn{x) is the nth Bessel function of the first kind. 
To proceed further, we consider three representative val- 
ues of s where the integration over momentum is conver- 
gent independent of the cutoff and we can set A — ;> c». 



g(o|ji'k-aki'k) _ g^lal.Vk,akvl] ^al^Vk g-Oki'k 



which results in 



(All) 



exp 



tXj2Fr.{A)at 



k^^O 



. _ -ji. Ti ' ri ' 1 



TT g-fki'k : g^k^k-ak-Uk 



k^^O 



7Ek=i0''k'"k 



g 2 Z^k^ 



X : exp 



■Y,F^^{A)if.l2) 



1. Subohmic case 

For s = —1/2 we can write 

Gi^HA) = / -Jo {k\r - s|) [1 - cos {kvA)] . 

(B3) 
Then, Eq. (I?!]) can be obtained by using the integral |29| 

— Jo(/3a;) [1 — cos(ax)] 



/o 



\fp^^~~o? + a arcsin ( ^ ) , a < /?, 



2 ' 



a>/3. 



(B4) 



where : (...) : denotes normal ordering. We can rewrite 
the argument of the first exponential in Eq. (|A12p since 



E-^4--tE E 



5r.,k5r,,k I <cr^ 



k^^O 



ij \k7^0 



ECl-l^X'^J, (A13) 



^,3 



where tj^.r (A) is defined in Eq. (IT5)) . Combining Eqs. 
(|M)) . (|A12I) . and (|A13p . we arrive at Eq. (fTS]) . 



Appendix B: The correlator t/rs (A) 



The correlator in Eq. ([T9| can be evaluated in the 
following way. Inserting Eq. (PH)) into Eq. p^ and 
using Eq. (|A1|) . we find 



5(«)(A) ^ (^/?jr'' V |k|2-2g.k.(r-s) 



k/0 



g»|k|t,A „ -1^ 



(Bl) 



2. Ohmic case 



For s = we can write 
,(i.)^.^ 1 rdk 



TTUJ^ Jq k 



Jo(fc|r-s|)[l-cos(fcuA)] 



(B5 



Q'r^K^) 



Then, Eq. ([23l) can be obtained by using the integral \^ 

f°° dx fa\ 

/ — JoiPx) [1 — cos(ax)] — arcosh — 6(a — P). 
Jo X \PJ 

(B6) 

3. Superohmic case 

For s = 1/2 we can write 

gi^\A)^ ^ / d/cJo(fc|r-s|)[l-cos(fcfA)]. 

(B7 



Then, Eq. (f25|) can be obtained by using the integral [3lJ 

dx MM [1 - cos(ax)] = i - ^£==- (B8) 



Appendix C: The correlator C/^s (A) 
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The correlator in Eq. (|T8)) can be evaluated in the following way. Starting with Eq. ((T3)) . we have 



[/(r,ti),/(s,i2)]--2*- 



LD 



^|k|2«sin[k.(r-s) + |k|«(ti-t2)], 



k#0 



which allows us to write 



- {[/(r, ii), /(s, t2)] + [/(s, ti), /(r, i2)]} = -2i 



Considering D = 2, wc have 



■■ ("/"^IJ^ 5^ |kr cos[k . (r - s)] sin[|kb(ti - h)]. 



kT^O 



(CI) 



(C2) 



U[f{r,h)J{s,t2)] + [f{s,h)J{r,t2)]} = --(—] f dkk^'+'Mk\r-s\)sin[kv{t,~t2)]. (C3) 

In Eq. (jCSp . we introduced an ultraviolet momentum cutoff A. To proceed further, we need to specify s. Below, we 
consider three representative values. 



1. Subohmic case 



For s = —1/2, the integral in Eq. (ICSP converges. Using the integral [32| 

9{\a\-\l3\) 



we have 



dx Jo (/3a;) sin(aa;) = sign(a) , 

v<^ ~ P 



i {[/(r,ii),/(s,t2)] + [f{s,t,)J{r,t,)]} = -'- (^) sign(ii - t,) ^^^LJ^=^=^. 
2 TT V^o/ ^u2|ti -i2p - |r-s|2 



Carrying out the time-ordered integration over ti and t2, we obtain Eq. 



(C4) 



(C5) 



2. Ohmic case 



For s = 0, we notice that 

-{[/(r,ii),/(s,i2)] + [/(s,ti),/(r,t2)]} = — - — — / -Jo(fc|r-s|)sinMti-i2)] 
2 TT V^o/ dt2 dti Jq k 

Using the integral [3j| 



— Jq(j3x) sin(ax) = sign(a) 



-0(|a|-|/3|) + arcsin^ ^(l/?l - |a|) 



we obtain 



i{[/(r,t0,/(s,i2)] + [/(s,ti),/(r,t2)]} = -- {-) ^-^ 
2 TT XuiQ J dt2 dti 



■e{\h-t2\-\v-s\) 



arcsmi ^^^ ) 0(|r - s| - |ti - ^2!) 
r — s 



(C6) 



(C7) 



(C8) 



Carrying out the time-ordered integration in ti and t2, we obtain Eq. (j24p . 
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3. Superohmic case 

Similarly to the olimic case, for s = 1/2 we write 



1 {[f(r,t,)J{s,t2)] + [fis,h)J{r,h)]} = --(-) 4r4- f dkJoik\r-s\)sm[kv{h~h) 

2 TT \oJoJ at2 ati Jq 



Using the integral in Eq. (IC4p . we obtain 

i{[/(r,ii),/(s,i2)] + [/(s,ti),/(r,i2)]} = --f-) ^^ 
2 TT Vwo/ dt2 dti 



9{\h-t2\-\v~s\) 
v/^2|i,_i2|2_|r-s|2 



Carrying out the integration the time-ordered integrations on ti and t2, we arrive at Eq. (1261) . 

I 



Appendix D: Single-qubit flipping probability 



since gill^ (A) = 0. 

Consider now the change of basis 



1 



1 



We can obtain a compact expression for the single- 
qubit-bath evolution operator in Eq. ([55)1 by following 
essentially the same steps show in Appendix [Xj The only 
formal difference is that summations over all qubits in 
the lattice have to be replaced by a term corresponding 
to a single qubit j. Thus, considering Eq. dH]), the result ^ff\y{ch allows one to write 
is 



t,> = ^ (I t,)x + 1 ;,).) , 



U,> = -^(|t,>x-U,).): 



1 



where 



Uj{^) = Xi :exp 



Xj = exp 



-y^r,(A)aJ 



(Dl) 



(t, |C/,(A)Uj)--iexp 



exp 



f-^^.(A)l 


;f^.(A) 


} 



(C9) 



(CIO) 



(D3) 
(D4) 



(D5) 



om^) 



since : e '^^•^i^^'^ -.^^^^e '^ ^r^i'^) [gee Eqs. (|A12I) and 
(D2) (Pl3l) ]. We can then write 



(t, |C/](A)| ;,)(;, |C/,(A)| t,) = \{2- exp [zXFr^iA)] - exp [-ai^r,(A)] } 

- i {2 - X,': exp [tX F,^ (A)] : -x ■ : exp [-^A F,^ (A)] :} , 



(D6) 



I 

which yields Eq. (j59p when the expectation value over the bath vacuum is taken. Now we can insert Eq. (|D1 

into Eq. ^ to obtain Eq. (l59l) . 
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